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An  algorithm  is  presented  for  the  rapid  evaluation  of  expressions  of  the  form 


m 


j= i 


at  multiple  points  xi,X2, ■  ,.rn.  In  order  to  evaluate  the  above  sum  at  n  points,  the  algorithm 

requires  order  0(n  +  m )  operations,  and  a  simple  modification  of  the  scheme  provides  an  order 
O(n)  procedure  for  the  evaluation  of  an  order  n  polynomial  at  n  arbitrary  real  points.  The  algorithm 
is  numerically  stable,  and  its  practical  usefulness  is  demonstrated  by  numerical  examples. 
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1.  Introduction 

In  this  paper,  we  present  an  algorithm  for  the  rapid  evaluation  of  expressions  of  the  form 

rn 

f1) 

j= i 

where  x  >  0,  a  =  jet  i .  ag,  •  • .  a,,, }.  J  =  { ,  3i< '  •  dm  }  are  two  finite  sequences  of  real 
numbers  and  >  0  for  all  1  <  /  <  in.  To  evaluate  the  sum  (1)  at  n  arbitrary  points  on  the 
real  axis,  the  algorithm  requires  a  number  of  arithmetic  operations  proportional  to 

(n  +  in  /f >(/•_.(  -))  (/o#yj (  —  ))".  (2) 

t  t 

where  e  is  the  precision  with  which  the  calculations  are  being  performed,  and  in  most  cases 
likely  to  be  encountered  in  practice,  the  estimate  (2)  can  be  reduced  to 

(n  +  in)  loiji(-)  (3) 

t 

(see  Observations  7.1.  7  2  below). 

The  evaluation  of  expressions  of  the  form  (i)  is  closely  related  to  several  classical  problems 
in  the  theory  of  computation.  For  example,  the  problem  of  rapidly  evaluating  a  polynomial 

P(/)  =  f \Pr0  (4) 

J  =  i 

at  in  different  points  is  readily  reduced  to  the  form  (1)  by  the  obvious  substitution  x  = 
lotj(l).  The  classical  algorithm  for  evaluating  (4)  at  in  points  has  an  asymptotic  complexity 
0{rn  ■  log2(rn))  (see,  for  example,  [1.2]),  making  (3)  a  moderate  improvement  over  previously 
available  results,  so  far  as  the  asymptotic  CPU  time  estimate  is  concerned.  On  the  other  hand, 
the  algorithm  of  the  present  paper  is  numerically  stable,  and  our  numerical  experiments  (see 
Section  8)  indicate  that  in  practical  calculations,  it  is  extremely  efficient,  making  it  a  method 
of  choice  whenever  expressions  of  the  form  (1)  have  to  be  evaluated  at  large  numbers  of  points. 

Remark  1.1.  Classical  algorithms  for  the  rapid  manipulation  of  polynomials  are  purely 
algebraic,  and  are  applicable  to  polynomials  over  a  wide  class  of  fields.  On  the  other  hand,  the 
algorithm  presented  here  is  based  on  approximation  theory  (i.e  it  relies  on  certain  facts  from 
real  analysis)  and  is  restricted  to  polynomials  over  the  field  of  real  numbers.  While  it  can  be 
generalized  to  certain  other  fields,  detailed  investigation  of  such  generalizations  is  outside  the 
scope  of  this  paper,  and  will  be  reported  at  a  later  elate. 

2.  Relevant  Facts  From  Approximation  Theory 

Suppose  that  n,b  are  a  pair  of  real  numbers  such  that  <1  <  h,  and  that  k  >  2  is  an  integer, 
(’hebvehev  node's  U  ./•„>,  ,/*  on  the  interval  [ei./i]  are  de'fineel  by  the  formula 


a  4-  b  n  —  b 


2  ;  +  1  rr 


For  a  function  /  :  [a.b]  — -  Rl .  we  will  denote  bv  Pjf  h  y  the  order  k  —  1  Chebyrliev  approxi¬ 
mation  to  the  function  /  on  the  interval  [«.//],  i.e  the  (unique)  polynomial  of  order  k  -  1  such 
that  P*k  j(ti)  =  /(/,)  for  all  i  =  1.2,  ,  A\  There  exist  several  expressions  for  the  polynomial 

P*h  j-,  and  the  one  we  will  use  in  this  paper  is 

*£»./(*)  =  £«/  (')-f(tj)  (C) 

]=  1 


with 


iij(t) 


nf=i 

n  ii^tj  -  t.y 


The  following  well-known  lemma  provides  an  error  estimate  for  Chebychev  approximations. 
It  is  the  principal  analytical  tool  of  this  paper,  and  can  lie  found,  in  a  somewhat  different  form, 
in  [3). 

Lemma  2.1.  If  f  €  e* [«.£»]  (i.e.  f  has  k  continous  derivatives  on  the  interval  [n.hj).  then 
for  any  t  €  [ft.  6), 


I  Pa.b.f(t)  ~  f{l)  |<  ^ 


( b  —  (/)*■' 
41 


(3) 


with 


M  =  max  |  /(*>{0  I  . 
te[a.h\ 


(0) 


Furthermore,  for  any  k  >  2  and  /  €  (a. 6], 


£“>(')  <2. 
;  =  i 


(10) 


ami 


(11) 


£  i  «>(0  i  ^  2  +  -  io,,{k). 

1= 1  ' 

In  the  present  paper,  the  above  lemma  will  lie  used  in  the  special  case  where  0  <  a  <  b. 
and  f(t)  =  e~";  *,  with  *■  >  0.  Under  these  conditions,  the  expression  (S)  assumes  the  form 


-  W)  l<  f; 


(12) 


and  the  following  lemma  provides  a  form  of  the  estimate  ( 12)  independent  of  7. 

Lemma  2.2.  If  under  the  conditions  of  Lemma  2.1.  f(t)  —  e~~,  l.b  =  2a.  and  a  >  0,  then 


rhrin  -fin  i<  tt 


(13) 


for  all  k  >  2  and  t  €  [«,/>]. 

Proof.  Obvously.  for  t  £  [a, 2d],  the  estimate  (12)  can  be  rewritten  in  the  form 


•/*  ak 


Differentiating  the  latter  expression  with  respect  to  7,  we  find  that  its  maximum  is  achieved  at. 


Now.  substituting  (15)  into  (14)  and  using  Stirling's  formula,  we  obtain 


P,k;,.f(t)~W)\<±  ■(£)*  -JF  ■«-*■"<£  ~  e-k  < 


3.  Exact  Statement  of  the  Problem 

In  the  description  of  the  algorithm  below,  we  will  assume  that: 

a)  d  =  {o  1 .  a>.  •  .n,„  ).  3  —  { J| ,  }.  x  =  •■,./„]  are  three  finite  sequences  of 

real  numbers. 

b)  The  sequences  3  and  i  are  monotonicallv  increasing. 

c)  Ji  >  0. 

d)  jq  >  0. 

e)  We  would  like  to  evaluate  tlie  suns 

m 

s„.A*k)  =  ]T'>.  • ,r 3‘  U7) 

i= 1 

for  all  k  =  1.2.  .  n  with  a  relative  accuracy  e  >  0,  i.e.  we  would  like  to  find  a  number 

)  such  that 

!  iv.f Uk)  -  rk.)  I 


n;=,  I'D  t 


for  each  k  €  (l.  «]. 


Remark  3.1.  As  h  as  been  mentioned  in  the  Introduction,  the  problem  of  evaluating  a 
polynomial  of  order  m  at  n  points  is  easily  reduced  to  the  form  (17).  Indeed,  suppe  e  that  a 
polynomial  of  the  form  (4)  has  to  be  evaluated  at  a  monotonicallv  increasing  finite  sequence 
4  points  ,l„.  It  can  be  assumed  without  a  loss  of  generality  that  0  <  <  1  for  all 

k  =  1.2.  .  n.  and  ne  will  introduce  a  new  vari  tide  .r  =  -loy(t).  and  denote  )  by  .1  >. 

Thus,  evaluation  of  the  polynomial  (4)  at  a  monotonicallv  increasing  finite  sequence  of  points 
has  been  reduced  m  evaluating  the  expression 
rn 

Yr,  (10) 


3 


at  a  monotonically  decreasing  finite  sequence  of  points  .r  =  |  x i .  Xj .  •  •  ■  .  x„  ] .  Finally,  by  revers¬ 
ing  tho  order  of  rh<>  ■'<■< i we  reduce  the  evaluation  of  the  polynomial  (4)  at  the  point* 
tj.tj.  to  the  standard  problem  fonmilateil  above. 

4.  Notation 

In  this  section,  we  will  introduce  several  definitions  to  be  used  in  t lie  description  of  the 
algorithm  in  Sections  5.  G  below.  Throughout  this  section,  we  will  assume  that  we  are  dealing 
with  the  problem  described  in  Section  3.  and  that  y  is  an  integer  whose  particular  value  is  to 
he  determined  later. 

We  will  denote  hv  M  the  smallest  integer  number  such  that 

.(■„  _ 


We  will  define  a  finite  sequence  [fq  ] .  /  =  1.2.  .  M  of  intervals  on  the  real  axis  bv  the  formula 

r'/  =  t'^'r  (-1 

Similarly,  we  will  define  a  finite  sequence  jV,}.  i  =  1,2,  ,  M  of  intervals  by  the  formulae 

A-1</<.U-1. 


\i  =  <>• 


For  any  /  =  1.2.  .  M.  we  will  denote  by  ,i,  the  subset  of  )  consisting  of  all  points  such 
that  it  €  l’,.  and  for  any  /  =  1,2.  .  M.  we  will  denote  by  x}  the  subset,  of  x  consisting  of  all 
points  ii  such  that  x}  €  V,. 

For  e.ic  h  i  -  1.2.  .M.m,  will  denote  the  number  of  elements  in  i,.  Similarly,  for  each 

/  =  1.2.  ..\f.ll,  will  denote  the  number  of  elements  in 

Remark  4.1.  Obviously.  < !< •] > < ’n<  1  in g  on  the  distributions  of  the  points  and  .r, .  the  .1/ 
i  an  be  fairlv  large.  However,  the  total  number  .1/  of  *uch  /  that  m,  ^  (l  is  bounded  by  m.  and 
the  total  number  .V  of -urh  /  that  u,  ^  (I  is  hounded  by  u.  For  obvious  reasons,  we  will  refer 
a*  einprv  to  interval'  fq ,  1 ’  sinli  that  iii,  =  0  and  iij  =  0  In  the  opposite  case,  the  intervals 
will  be  referred  to  as  non-empty. 

For  each  i  =  1.2.  .  M.  and  /  —  1.2.  .  </.  we  will  denote  by  V*  the  j-fh  (Hiebychev  node  on 

the  interval  Ut. 

Similarly.  for  each  /  =  1.2.  .  M.  and  j  =  1.2.  .</.  we  will  denote  by  x'  the  j-r  h  <  Mteby- 

f  h e v  node  on  the  interval  1 " , 

For  each  k  =  1.2.  .  M.  and  >  'Ut  h  that  j,  c  f’k.  we  will  define  the  finite  sequence  \uk  \.  j  = 

1.  2.  .  i  bv  r  lie  f.  .rmn’.a 


'1  -I 


For  each  k  =  1. 2.  M.  and  i  —  1.2,--,  we  will  define  a  real  number  u*  by  the  formula 

53  E  °j  "?.>•  (2 

■t;  6f't 

Observation  4.1.  Due  to  Lemma  2.2.  the  expression 


>.*(o  =  E  "f., 


can  be  viewed  as  an  approximation  to  the  function  e  Furthermore,  for  any  t  €  ;0.  oc]. 

I  of'(0  -e-1'1  |<  1.  (■ 

4  ' 

Combining  (24).  (25).  (26)  with  the  triangle  inequality,  we  easily  see  that  the  sum 

».«>  -  t <  '»£  E  «>• <■ 

1  =  1  i  =  l  3,£l'k 


=  E  °>  E'E 

.y,ef*  .=i 

can  be  viewed  as  an  approximation  to 

E 


and  that 


:  OfcIO  -  21  ■  «  j‘  '  1$  4  E  I  rL  I  ■ 


Furthermore,  combining  (11)  with  (24)  and  using  the  triangle  inequality,  we  obtain 


E 1  "f E  fi  nJ  i  Ei  "o  i)  <  (‘2  +  5  E  I  !  •  (: 

i=l  j€l-\  ,=  l  i,€l’k 


<  liven  I:  =  1.2.  .  M.  and  i  =  1.2. 


.7 .  we  will  define  a  real  number  /*  by  the  expression 


(42) 


Observation  4.3.  It  is  easy  to  see  that  if  x  €  U,  and  3  €  Yj  with  j  <  then 


r -3 


<  e. 


Similarly,  if  .r  €  U,  and  3  6  \  'j  with  j  >  //,.  then 

I  e',  J  -  1  |<  e.  (43) 

Fuerthermore.  for  any  i  —  1.2 .  .M  —  1. 


<  2  /y</)(-).  (44) 

e 

In  other  words,  given  .r  £  U,  and  3  €  I one  of  three  possible  situations  obtains: 

a)  ./  <  i-',-.  In  this  case.  e~J‘*  can  be  approximated  by  0  with  a  precision  > . 

b)  j  >  //,.  In  this  case.  e~  '*  can  be  approximated  by  1  with  a  precision  t. 

cl  i'i  <  j  <  /i,.  In  this  case,  e~  1  *  can  not.  be  approximated  by  either  0  or  1.  However,  the  total 
number  of  indices  j  for  which  this  situation  obtains  is  bounded  by  2  loj2{  f ).  independently  of 
x.  3.  or  i . 

5.  Informal  Description  of  the  Algorithm. 


We  will  illustrate  the  idea  of  the  algorithm  on  a  Mtnplified  example.  Namely,  we  will  assume 
that  3t  £  U\ ,  i.e. 


~  <  i  <  ? 

9  b  ^ 


for  all  i  =1.2,  ,  rn .  and  .r }  6  \\*  i.e. 


o 


<  <  -r 


(45) 


(4G) 


for  all  j  =  1. 2.  •  .  n. 

Consider  the  function  n~i  £  with  3  G  U\ .  x  &  \  \  .  Fixing  x  and  viewing  e~  J  J'  as  a  function 
of  3.  we  construct  its  (/-point  Chebychev  approximation  c'l(3)  on  the  interval  U \.  Due  to  (G). 


'■  '1(3)  =  'M  ^  '  J’  *  • 
s  =  1 


(47) 


with  the  functions  iij  debited  bv  (7).  and  the  coefficients  3‘  defined  in  Section  4.  According 
to  Lemma  2.2. 


/uul,  giv^ii  a  fixed  precision  t.  we  can  choose  >\  -  2  / # **/ 1 1  f).  md  in  a]j  ^uh^/pi^nr  r;t  j'nLiUoj].- 
replace  e~  *  r  with  ('oiiihiuim;  (4$)  with  the  triangle  ineipialitv,  we  obtain  the  ■■-rimaTi 


V  ,v . .  _  V' 


i 

17  ■  £  !  1 


for  any  ./•  g  (0.  ~-ccl.  und  *  1  <  1  e  to  (-45).  the  latter  can  l>e  rewritten  in  the  form 

<f  rn  .  in 

w,.,  £„r, ,-*'!<  ±  £!„,(.  tm 

,=1  J  —  1  J  =  1 

with  the  coefficients  '.'t  •  t  ,/  defined  by  the  formula 

m 

<■',  =  Yr»r  (51) 

;  =  i 

Now.  instead  of  evaluating  (17)  at  each  of  t he  points  .r,.  we  start  with  evaluating  the 
coefficient*  =  1.2.  •  .</.  which  is.  obviously,  an  order  ()(m  7)  procedure.  After  that,  we 

evaluate  the  expression 

t‘.,  (52) 

i=l 

for  all  /.'  =  1.2.  .  n.  which  i>  an  order  Oln  7)  ]>rocedure  (evaluating  a  7- term  expansion  at  n 

points) .  Thus,  the  total  operation  count  becomes  (J{  (  n  ~  m )  ■  7 ) .  Due  to  ( 4 J ) .  in  order  to  obtain 
a  relative  accuracy  , .  7  has  to  be  of  the  order  ).  and  we  have  reduced  rite  computational 

complexity  of  evaluating  (17)  from  ()(n  m)  to 

0((n  +  m)  /07.1(h).  (53) 

An  alternative  approach  would  be  to  calculate  the  coefficients  v,  for  i  =  1.2.  ■■.7  (order 
m  ■  7  operations),  evaluate  the  expression 

TV,  (54) 


for  all  /.•  =  1.2.  . -/  (order  7*  operations  ).  and  interpolate  the  expression  (17)  from  the 

<  hebyi  lev  node-  .1  | . ,/  l.  .  .r,1.  t  o  t  lie  points  . r  | .  .cj.  •  ■../•„  (oiaier  n  7  operations).  The  result¬ 
ing;  CPI’  time  estimate  in  this  case  is 

Ol  (n  -  m  )  7)  a-  TM7-)  =  Of  (/i  ~  m)  /o,/.,  (  -  )  +  O  ( ( />>>/ ,  (  -  ) ) J ) .  (55) 

'.  ‘  ( 

which  is  not  substantially  different  from  (53). 

When  the  points  d|..+2.--  .  i,„  and  .i|..r>.  do  not  satisfy  the  inequalities  (45). 

(4Gj.  the  above  approach  -  an  not  be  used  in  'U<  h  straightforward  manner.  However.  for  an\ 
i.  j  €  1 ,  .\ /  .  Lemma  2.2  can  be  used  separately  on  each  of  the  intervals  with  the  results 

'  ombined  to  of, tain  an  approximation  to  (17).  Thi'  i-  done  in  the  following  section,  resulting 


(5G) 


0{[  n  -r  in  )  /(.';/(  —  )  —  ()(n  (/('</(  -  ) )') 

<T  /' 

algorithm  tor  evaluating  (17)  at  it  points  with  a  relative  precision  t. 

6.  Detailed  Description  of  the  Algorithm 

Algorithm 
Stage  1. 

Comment  Choose  parameters  and  perform  geometrical  preprocessing. 

<  'house  precision  ,  to  he  achieved.  Set  •/  =  2  /<•«/(  i).  ( instruct  the  intervals  fft ,\' 
sets  V, .  with  /  =  1. 2.  ■  .  .1/. 

Stage  2. 

Comment  On  each  of  tie'  non-empty  intervals  Uk.  evaluate  the  coefficients  uk  in  the 
'ions  (27).' 

Step  1. 

Comment  'Set  all  coefficients  m*  to  zero.] 

do  k  =  1 .  A /  —  I .  C  0 
do  i  =  1. 1/ 

set  nk  to  zero. 

end  do 
end  do 


Step  2. 

Comment  For  •  n  h 

tie.  »»*.; 

do  k  =  I .  A /  —  1 .  h  vt  ,1 
do  i  =  1.7 

do  -j  e  t\- 

Evaluate  ./  Via  t-rmuia  i  A3  i  and  add  the  l>l"du<t  .  </*  r,>  ,/ 

y. '  1  1  j  <  t 

end  do 
end  do 
end  do 


h  of  rhe  |]on-ciii|.ty  intervals  JTk.  evaluate  n, 


Stage  3. 


and  the 


>*.n  pan- 


id  it  to 


0 


Comment  Evaluate  /*  via  formula  (31)  for  ail  !:  -  1.2.  .  3/  -uri.  ri.  ,» 


do  =  1.3/  —  1 .  .<-fr  =  0 
do  i  =  1 .  ij 

evaluate  the  expression  /*'  =  — '/=i  "i  1;  ‘ 

end  do 
end  do 


Stage  4. 

Comment  r«U  each  j  =  1.2.  ,  n.  evaluate  fj  via  formula  (32). 

do  I:  ~  1.3/-  l..rk  3=  >3 
do  ./  tl'( 

evaluate  the  expression  /,  =  2^/=1  r,1'  1.1*' 

end  do 
end  do 


Stage  5. 

Comment  Fair  earli  /.■  =  1.2.  .3/  and  each  a  ,  6  1 .  use  observation  4.3  to  walnut 

sum  Sr  =  3Z,(|c  ny  o }  ■  Add  the  result  to  /,.  eonduding  the  calculat  ion.! 

Stej)  1 

Comment  'Evaluate  St.; 

<,‘t  r, 


Step  2 

Comment  Evaluate  St  reeursively  foi  /,-  =  2.3.  .3/.] 

do  /•  =  2.  3/.  ^  0 

evaluate  Sr  via  the  tonnnla  Sr  =  S/._|  -j-  ^  <  a\  1 

end  do 


Step  3 

Comment  F>r  til  k  -  1.2.  .3/.  and  all  /  siidi  that  r,  e  \  \  .  a. Id  Sr  to  ,i , .  . . uu  lud in 

e  i  a  i  ,,i  r  a  it . 


ID 


if}  iff 


do  x,  €  \  \ 
add  5;.  to  /,. 

end  do 
end  do 


7.  Complexity  analysis 

Stage  Operation  Explanation 

number  count 

St, ice  1  Ola  -+-  m)  Each  of  the  points  .fi.do.-  ■ .  ,f((,  is  assigned  to  a  single 

interval  Ut.  Each  of  the  points  .r(  ..to.  •  •  -  .  .t'„  is  assigned 
to  a  single  interval  1 '■ . 

race  2 

tep  1  0(M  ■  q)  Each  of  the  coefficients  u*'.  with  k  =  1.2.  ■  .M. 

and  /  =  1.2.  •  •  ,  </  is  set  to  zero. 

Step  2  G(m  i]2)  Each  of  the  points  d| .  ,d),  ■  . 

contributes  to  the  coefficients  u1'  t  with 
j  =  1.2,  -  •  •  .</,  and  evaluating  each  of  the  coefficients 
«*,  requires  order  q  work  (see  (23))  . 

Stage  3  (){M  q2  (k))  The  sum  (31)  has  to  be  evaluated 

at  </  nodes  j-f.j-*.  ■  •  ■  x*  on  each  of  non-empty 

intervals  1  '| .  1  2.  • .  1  ,\t. 

and  on  the  I:  -  Ih  interval,  it  contains 

/'*•  —  f-'r-  terms.  However,  due  to  (4-1)).  p*.  —  i/j.  <  t 

for  all  k  =  1,2.  ,.V. 

Stage  4  G[n  >/*)  The  expression  (32)  has  to  be  evaluated  lot- 

each  of  the  points  ,j-|  .  ,r  >.  ■  ■  • .  .r„.  and  evaluating 

each  of  the  coefficients  r*  requires  order  </  work  (see  (33 

Stage  h 

Step  1  Ofm)  The  sum  5i  =  V)  )  gr,  o,  contains  no  more  than  m  terms 

■-r,‘P  2  0[n  -+-  m)  The  total  number  of  non-empty  intervals  1'  In  bounded  b\ 

and  the  total  number  of  coefficients  n) 
is  bounded  bv  in. 


Stop  3 


(){u\ 


Each  of  the  numbers  f}  is  amended  once. 


Summing  up  r h . ■  ('PE  rimes  for  all  stages  above,  we  obtain  the  following  time  estimate: 

T,..t,a  —  ■'  m  ~  a  n  e  m  ij~  •  I  n  M  </'  I  *  "j  (  -  ) .  (57 ) 

where  rin>  coefficients  t.n. depend  on  the  computer  system,  language,  implementation, 
etc.  However .  M  <  m.  .uni  /  f  i.  ami  the  estimate  (57)  assumes  the  form 

r,„Mi  =  a  hi  i  I  '  - /•  n  (/",/(-  (58) 

r  r 

The  estimate  |5S))  i'  iin  lep'Ui1  lent  e.f  the  locations  of  the  points  in  /?'.  ami  <loes  not 

depend  on  any  pie.  imputed  -lata.  The  following  two  observations  re<luce  it  to 

T,  ,/.,t  =  Phin  —  „  i  /o,/i  -  j  (59) 

for  many  problems  of  practical  interest. 

Observation  7.1.  The  term  h  m  ./-  in  (58)  is  associated  with  the  Stage  3  of  the  algorithm 


ami  t he  grossly  pessimistic  estimate 

,T7  <  M  <  in.  (GO) 

Aer-or-ling  to  j2i>i. 

M  -  1  <  - )  =  /'Vil in )  f  /“:/:(•*' i ■ )  +  (  -  )•  (01) 

*  t 

Normally,  when  calculations  are  performe<I  on  i  physical  computer,  the  exponential  in  the 
binary  representation  of  a  real  number  is  bounded.  ami  we  will  <lem>te  this  liouml  by  A.  It 
immediately  follows  that  in  all  cases. 

Tl  <  M  <  3  Inn  A  A)  -  1.  (02) 

ami  the  estimate  i57)  beomi.s 

T,„i„t  =  n  ni  —  h  it  <  m  -h  , I  n  ’i~  +  •  l<"i(A)  i('  (03 ) 


Observation  7.2.  The  terms  r  m  q-  and  d  n  q~  in  (57)  tie  associate!  with  tlie  Stages 

2  ami  -4  of  tin*  algourhni.  ami  with  rh-*  fact  that  in  order  to  evaluate  each  .4  the  eoetficients 

a ^  l  (or  /  *  (  i.  a  q  -  1-r.um  pro* bu  r  of  rlie  form  (24)  (or  (33))  has  to  be  evaluated.  <  Ibvionslv. 
tin*  eoefhid'-nfs  uk  depend  only  on  tin*  distubufion  of  rlie  points  t,.  and  n-a  n  that  of  ./•,  or 
a,.  Similarly,  the  roefticients  rk  depend  only  on  tin*  distribution  . .f  the  points  and  not  on 

that  ..f  it  or  ■»,.  Therefor**,  for  fixed  distributions  of  i,  and  r,.  the  >  0,4*1.  ants  .•*  (  can  be 

pr*-i  ompur**d  and  stored.  r*-dueiiig  the  total  ('PI-  rime  estimate  to 

12 


Ttoiai  ^  ,f  fli  7  +  n  7  -t-  +<*  •  It'il(A)  7“  •  /mi/(  -).  (04) 

t 

However,  q  s.  and  U»j(A)  is  fixed  for  given  computer  system  and  language.  Thus,  when 


Ttol„i  (n  ■  in  +  b  n)  ■  q 


8.  Numerical  Results 

A  computer  program  has  been  written  implementing  the  algorithm  of  this  paper.  The 
calculation  is  performed  in  two  stages,  each  implemented  by  a  separate  subroutine.  During 
the  initialization  stage,  the  coefficients  »/*  »•*'  are  ev.aluate<l  for  given  distributions  of  points 

■  fj.  .  j’i  ..»■.>.•■  .  .r„  (see  Observation  7.2).  During  the  second  stage,  t  lie  sums  (17) 

are  evaluated  for  a  given  set  of  weights  oj.m.j.  .  o,„. 

Remark  8.1.  It  is  clear  from  Tables  1.  3.  a  that  the  hist  stage  (initialization)  tends  to 
be  several  times  more  expensive  than  the  second  (evaluation).  However,  in  most  applications 
the  algorithm  has  to  be  initialized  once,  with  subsequent  repeated  evaluation  of  the  sums  (17) 
for  varying  sets  of  weights  .  o,„.  This  situation  is  similar  to  that  encountered  for  the 

Fast  Fourier  Transformation. 

The  program  has  been  applied  to  a  variety  of  situations,  and  three  such  examples  are 
presented  in  this  section,  with  the  computations  performed  on  a  VAX-SGOO  computer.  In  each 
case,  we  performed  the  calculations  in  three  way>:  via  the  algorithm  of  the  present,  paper 
in  single  precision  arithmetic,  directly  in  single  precision  arithmetic,  and  directly  in  double 
precision  arithmetic.  The  first  two  calculations  were  used  to  compare  the  speed  and  precision 
of  the  algorithm  with  that  of  the  direct  calculation.  The  direct  evaluation  of  the  field  in  double 
precision  was  used  as  a  standard  for  comparing  the  accuracies  of  the  first  two  calculations.  In 
all  cases,  we  set  c  =  1()-8.  and 

m  =  n  =  10  ■  2fc.  (GO) 

with  I:  varying  from  1  to  S. 

Tables  1.  3.  5  contain  the  CPU  timings  for  the  examples  1.  2.  3  respectively.  Following  is 
a  detailed  description  of  the  entries  in  these  Tables. 

n  -  the  number  of  points  at  which  the  sum  ( 1 )  is  being  evaluated. 

T„„,  -  tin-  initialization  time  of  the  algorithm. 

T.,1,1  -  the  CPU  time  required  by  the  algorithm  once  it  has  been  initialized. 

T,Ur  -  the  CPU  time  required  by  the  direct  calculation. 

Tables  2.  4.  G  contain  the  accuracies  for  the  examples  1.  2.  3  respectively.  In  the  description 
of  the  entries  of  these  tables  below.  5*  denotes  the  sum  (1)  at  the  point  .iq.  as  evaluated  directlv 
in  double  precision.  S'l'r  denotes  the  sum  (1)  at  the  point  ,/q.  as  evaluated  directlv  in  single 
pp'ri'ion.  and  denotes  the  '•urn  (1)  -it  the  point  .iq.  as  evaluated  in  single  precision  via  the 


algorithm  of  the  present  paper.  Following  is  a  detailed  description  of  the  entries  in  the  Tables 
2.  4.  G. 

n  -  the  number  of  points  at  which  the  sum  (1)  is  being  evaluated. 

<V""r  -  the  maximum  error  produced  by  the  algorithm  at  any  point.  It  is  defined  by  the  formula 


<  n mix 

,d.j 


max 

i  <t<n 


(07) 


'  r(u>  maximum  error  produced  by  the  direct  calculation  at  any  point.  It  is  defined  by  the 
formula 


y:;:;'  =  max  i  s£ir  -  s, 

'  1  </.</! 


(GS) 


nni  i  .rel 


■  (  h  -  the  maximum  relative  error  produced  by  the  algorithm  at  any  point.  It  is  defined  bv 
the  formula 


<;  nmx.rel 


•*/«/ 


max 


*1*  -  5t 


l<k<n  Or 


(GO) 


b”' "J 1  n  -  the  maximum  relative  error  produced  by  the  direct  calculation  at  any  point.  It  is 
defined  bv  the  formula 


Sfr  -  5r 


~  'n<*X 


1  <*■</!  |  -5/. 


ro) 


-  rlie  relative  error  as  defined  in  Section  3  as  produced  by  r lie  algorithm.  It  is  given  by  the 
formula 


<rr.l  . — ■  t—  ]  I 

dal9  =  - 


£2=1 1  Sk 


-1) 


y)'  -  the  relative  error  as  defined  in  Section  3  as  produced  by  the  direct  calculation.  It  is  given 
bv  the  formula 


£2=,  \st-st 


E2=,  I  ^ 


Following  is  a  detailed  description  of  the  three  example- 


Example  1.  In  this  example,  the  points  d(.  i>.-  .  J„,  and  ,.r„  were  defined  lay 

the  formulae 


i, 


m  -  1 


-  1). 


a  = 


(/•■  -  1). 


73) 

<4) 


[4 


and  the  weights  »i,Q2,  •  ■  ,am  were  generated  randomly  on  the  interval  [0, 1].  Here,  by  "direct 
algorithm”  we  mean  a  straightforward  implementation  of  the  formula  (17).  The  results  of  this 
set  of  experiments  are  summarized  in  Tables  1,  2. 

Example  2.  In  this  example,  the  points  •  ■  ■  ,/?m  and  xi,x%,  ■  ■  ,xn  were  generated 

randomly  on  the  interval  [0,5],  and  the  weights  ai,a2,  ■  ,c*m  were  generated  randomly  on  the 
interval  [0.  l].  Again,  by  "direct  algorithm"  we  mean  a  straightforward  implementation  of  the 
formula  (17).  The  results  of  this  set  of  experiments  are  summarized  in  Tables  3,  4. 

Example  3.  Here,  we  evaluate  a  polynomial  of  order  n  at  a  collection  of  randomly  gener¬ 
ated  points  on  the  interval  [0, 1],  The  coefficients  of  the  n  —  th  order  polynomial  are  randomly 
distributed  on  the  interval  [0,  l].  In  this  example,  the  direct  evaluation  of  the  polynomials  is 
performed  via  the  Horner's  rule  (see,  for  example,  [3]),  and  the  algorithm  of  this  paper  is  ap¬ 
plied  via  the  formula  (1).  The  results  of  this  set  of  experiments  are  presented  in  Tables  5, 
6. 

The  following  observations  can  be  made  from  the  Tables  1-6.  and  are  in  agreemant  with 
the  results  of  our  more  extensive  experiments. 

1.  In  all  cases,  the  accuracy  produced  by  the  algorithm  of  the  present  paper  is  comparable  o 
that  obtained  by  the  direct  calculation.  For  large  n,  the  algorithm  tends  to  be  slightly  more 
accurate. 

2.  The  CPU  times  and  accuracies  produced  by  the  algorithm  are  virtually  independent  of  the 
distributions  of  points  a,-,/?i,xjt  in  R1 . 

3.  When  used  for  evaluating  expressions  of  the  form  (17),  the  algorithm  becomes  faster  than 
the  direct  calculation  at  n  —  m  <  20,  if  the  initialization  time  is  ignored.  If  we  include  the 
initialization  time,  the  break-even  point  is  between  n  =  m  —  40  and  n  =  m  =  60. 

4.  When  used  for  evaluating  polynomials,  the  algorithm  becomes  faster  than  the  direct  calcula¬ 
tion  at  roughly  n  =  m  =  40,  if  the  initialization  time  is  ignored.  If  we  include  the  initialization 
time,  the  break-even  point  is  roughly  n  =  m  =  300. 
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Table  1 

Example  1  :  Timings 

T,  ,„t 

Talg 

T,itr 

0.0112 

0.0015 

0.0081 

0.0369 

0.0042 

0.0318 

0.0802 

0.0092 

0  1278 

0.136 

0.0165 

0.5202 

0.218 

0.0283 

2.069 

0.333 

0.0468 

8.368 

0.484 

0.0784 

33  25 

0.727 

0.137 

133.58 

Table  2 

Example  1:  Accuracies 

tmax 

°a‘g 

tmaz 

°dir 

cmax,rcl  rmaz,rcl 

°alg  °dtr 

cre( 

°atg 

br,d 

dir 

412E-06 

.638E-06 

.383E-06  243E-06 

.359E-07 

.556E-07 

.  179E-05 

.219E-05 

.459E-06  .418E-06 

.783E-07 

.960E-07 

408E-05 

.688E-05 

.623E-06  .671E-06 

.  100E-06 

.  169E-06 

.  138E-04 

.262E-04 

.825E-06  .  1I6E-06 

.173E-06 

•326E-06 

.378E-04 

.873E-04 

.597E-06  195E-05 

.238E-06 

.550E-06 

922E-04 

.231E-03 

.  103E-05  .258E-05 

.279E-06 

.699E-06 

.273E-03 

.740E-03 

.841E-06  .473E-05 

.420E-06 

•  114E-05 

522E-03 

.233E-02 

.880E-06  886E-05 

.407E-06 

.181E-05 

o\ 


Table  3 

Example  2  :  Timings 

T,  ,i,t 

T„tg 

T,i,r 

0  0097 

0.0011 

0.0083 

0.0275 

0.0033 

0.0332 

0.07G8 

0.0089 

0.1328 

0.126 

0  0157 

0.536 

0.210 

0.0271 

2.12 

0.326 

0.0455 

8.50 

0.497 

0.0784 

34.12 

0.698 

0.1351 

138.34 

Table  4 

Example  2:  Accuracies 

k  max 
aig 

smaz 

°di  r 

*  max,  rtl  cmaxjtl 

°alg  °dir 

crel 

°alg 

°d,r 

312E-0G 

.  164E-06 

.  160E-06  .192E-06 

.268E-07 

.141E-07 

109E-05 

.270E-05 

.405E-06  .316E-06 

.501E-07 

.124E-06 

354E-05 

.364E-05 

.658E-06  .860E-06 

.832E-07 

.857E-07 

835E-05 

1G1E-04 

.8G0E-06  .845E-06 

.  lOOE-OG 

.194E-06 

108E-04 

298E-04 

.974E-0G  .  158E-06 

.115E-06 

.173E-06 

G06E-04 

128E-03 

.824E-06  .288E-05 

.185E-06 

.389E-06 

33GE-03 

657E-03 

.914E-06  423E-05 

.207E-06 

.100E-05 

451E-03 

149E-02 

.827E-06  .799E-05 

.348E-06 

.115E-05 

I 

Table  5 

f 

s" 

Example  3  :  Timings 

"1  n 

C" 

Tinit 

Talg 

Td,r 

i  20 

0.0135 

0.0015 

0.0013 

1  40 

0.0435 

0.0052 

0.0047 

s'  80 

0.0948 

0.0109 

0.0179 

t  160 

0.1327 

0.0172 

0.0729 

>'  320 

0.222 

0.0286 

0.2779 

^  640 

0.306 

0.0445 

1.101 

^  1280 

0.422 

0.0718 

4.54 

r  2560 

/ 

0.664 

0.1322 

18.37 

1 

Table  6 

Example  3:  Accuracies 

i 

r. 

cma* 

°alg 

f.maz 

°dir 

cmaz.rtl  <-max,rei 

°alg  ddtr 

Ktg 

it-ti 

°dir 

r.  20 

772E-06 

.260E-05 

.301E-06  .305E-06 

.673E-07 

.227E-06 

£  40 

218E-05 

.275E-05 

.396E-06  .337E-06 

.955E-07 

.121E-06 

80 

518E-05 

.554E-05 

.528E-06  .210E-06 

.127E-06 

.136E-06 

j  160 

615E-05 

.462E-05 

.671E-06  .304E-06 

.766E-07 

.576E-07 

320 

103E-04 

.232E-05 

.851E  06  .387E-06 

.646E-07 

.146E-06 

640 

224E-04 

.295E-04 

.832E-06  .422E-06 

.680E-07 

.894E-07 

1280 

615E-04 

.360E-uJ 

.745E-06  120E-05 

.949E-07 

.555E-06 

^  2560 

757E-04 

.120E-02 

.862E-06  . 191E-05 

590E-07 

.932E-06 

* 

>  ' 
t 

J 

J 

J 

J 

